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Abstract 

The present study applies the power spectral analysis technique to understand the diffusional 
anomaly in liquid silica, modeled using the BKS potential. Molecular dynamics simulations 
have been carried out to show that power spectrum of tagged particle potential energy 
of silica shows a regime with l/f a dependence on frequency / which is the characteristic 
signature of multiple time-scale behaviour in networks. As demonstrated earlier in the case 
of water (J. Chem. Phys., 122, 104507 (2005)), the variations in the mobility associated 
with the diffusional anomaly are mirrored in the scaling exponent a associated with this 
multiple time-scale behaviour. Our results indicate that in the anomalous regime, as the 
local tetrahedral order decreases with temperature or pressure, the coupling of local modes 
to network reorganisations increases and so does the diffusivity. This symmetry-dependence 
of the vibrational couplings is responsible for the connection between the structural and 
diffusional anomalies. 
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I. INTRODUCTION 



Tetrahedral network-forming liquids are characterised by strong, local anisotropic inter- 
actions which impose a local tetrahedral order with four nearest neighbours, in contrast 
to the local icosahedral order characteristic of the random closed-packing structures seen 
in simple liquids. The interactions in simple liquids are dominated by strong, short-range 
repulsions and weak isotropic attractions; typical examples are the rare gases and liquid 
metals^. The best known example of a tetrahedral liquid is water where each oxygen atom 
can form four tetrahedrally disposed hydrogen bonds with its neighbours. At room temper- 
ature, the hydrogen bond energy is of the order of 5-10/csT and the hydrogen bonds in water 
can break and reform on timescales of the order of picoseconds resulting in a fluctuating, 
random three-dimensional network of water molecules^ 3 . A second important example is 
liquid silica where the SiC>4 tetrahedra are linked to form a random network^. Interactions 
in Si02 have a mixed ionic-covalent character with Si-0 bond energy of approximately 450 
kJ mol -1 such that in the temperature regime above 2000K where the liquid phase is stable, 
it is of the order of 25k B T or less 6 . The liquid phases of Si, Ge, Sb, Bi and Ga are expected 
to be tetrahedral as are those of several ionic compounds such as BeF 2 and intermetallics 
such as InSb, GaAs, GaP, HgTe, CdTe and CdSei. 

The liquid state thermodynamic and kinetic properties of tetrahedral liquids are in many 
respects anomalous when compared with those of simple liquids 2,5 ' 8 . The best known of the 
thermodynamic anomalies is the density anomaly shown by the existence of a temperature of 
maximum density (TMD). Unlike simple liquids where density decreases monotonically with 
temperature, in tetrahedral liquids, increasing thermal kinetic energy over certain density 
regimes results in progressive destruction of the network, consequent compaction and a 
resulting increase in density. At atmospheric pressure, the TMD of water occurs at 279K 
whereas that of silica occurs at approximately 5000K. The behaviour of response functions of 
tetrahedral liquids, such as the isothermal compressibility (k) and the isobaric heat capacity 
(C p ), can also be anomalous. In simple liquids , both C p and n decrease monotonically 
with decreasing temperature but, for tetrahedral liquids, both k and C p as a function of 
temperature along an isobar can show minima. Thermodynamic arguments can be used 
to show that the existence of a density anomaly implies the existence of compressibility 
and heat capacity anomalies 2 . In addition to the thermodynamic anomalies, tetrahedral 
liquids also show kinetic anomalies. For example, the diffusional anomaly corresponds to a 
regime in which the diffusivity increases as a function of density, in contrast to simple liquids 
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where diffusivity decreases with density due to increasing steric hindrance. Qualitatively 
this anomalous behaviour of the mobility can be understood in terms of the disruption of the 
tetrahedral network by increasing pressure or density which facilitates translational motion 
of the water molecules. The phase diagram of tetrahedral liquids can show a number of 
unusual features, such as a negatively sloped melting line on the pressure-temperature phase 
diagram. There is also increasing evidence that tetrahedral liquids can show polyamorphism 
i.e. the existence of low and high-density liquid or glassy phases. 

Over the past decade, a better understanding of the connection between the various 
anomalous properties of tetrahedral liquids has been obtained from a combination of exper- 
imental and simulation studies on water, and to a lesser extent, silica. A key innovation 
has been the introduction of order metrics to quantify the nature as well as the extent of 
structural order present in such liquids*^. The extent of local tetrahedral order present 
around a tagged tetrahedral centre (e.g.O atom in water and Si atom in Si0 2 ) is gauged by 
a tetrahedral order parameter q defined as 



where ipjk is the angle between the bond vectors and where j and k label the four 
nearest neighbour atoms of the same type e.g. the qo order parameter in water is defined 
using the positions of the four nearest oxygen atoms&ii. The translational order parameter, 
r, was defined as 



where £ = rp 1 ^ 3 and £ c = 2.83^. The translational order parameter, as defined above, will 
capture the extent of pair correlations present in the liquid and can be thought of as a type 
of density-ordering, since r will increase as the random close-packing limit is approached. 
The structure and dynamics of tetrahedral liquids can be understood in terms of an inter- 
play between short-range tetrahedral ordering and translational or density ordering which 
promotes local icosahedral order seen in simple liquids. At a given temperature, q will show 
a maximum and r will show a minimum as a function of density; the loci of these extrema 
in the order define a structurally anomalous region in the density-temperature (pT) plane 
where the tetrahedral and translational order parameters are found to be strongly corre- 
lated. The region of the density anomaly, where dp/dT > 0, is bounded by the structurally 
anomalous region. The diffusionally anomalous region (dD /dp > 0) encloses the boundaries 
of the structurally anomalous region in the case of silica. In water, on the other hand, the 
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boundaries of the diffusionally anomalous region lie between those of the structural and 
density anomalies. Thus one can think of tetrahedral liquids as displaying a cascade of 
anomalies where structural, diffusional and density anomalies occur consecutively as the 
degree of tetrahedral order is increased. Our understanding of the molecular level processes 
which give rise to this cascade effect is still incomplete. An important step in this direc- 
tion is the demonstration of the cooperative origin of low-density domains in water which 
suggests that the probability density of centres of high local tetrahedral order must reach a 
certain threshold value before cooperative effects can precipitate the formation of ramified, 
low density clusters of water molecule a 11 ^ 3 . 

In a series of recent papers, we showed that power spectral analysis could be used to 
understand the connections between the local order, hydrogen bond network dynamics and 
the kinetic anomalies in wate r 14 i 15 i 16 i 17 i 18 i 19 . The power spectral density of an observable 
A(t) as a function of time t is defined as follows 20 : 



The power spectrum contains the same information as the time correlation function of the 
observable but it highlights its behaviour over a range of frequencies. The multiple time-scale 
dynamics of networks is reflected in the power spectrum, S(f), as a l// a dependenc e) 21 ^ 22 . 
In the case of water, where the network is formed by the coupling of individual molecules 
by hydrogen bonds, the observables which are sensitive to local motions as well as network 
re-organisations were found to be the local tetrahedral order parameter and the tagged 
molecule potential energy. Fluctuations in these observables generated power spectra with 
three key features: (a) a local vibrational peak corresponding to the librational modes in 
water; (b) a multiple time scale region with l/f a behaviour and (c) a crossover to white 
noise at lower frequencies identifying the time scales for decay of correlations in the network. 
The diff usivity was shown to be strongly correlated with the exponent a of the 1/ f a region 
of the power spectra of the tagged particle potential energy fluctuations. Thus the power 
spectrum has been shown to provide a direct dynamical measure in the scaling exponent a of 
the degree of coupling of the librational modes to the network vibrations which complements 
the static structural measures, such as the coordination number distributions, used in earlier 
studies. Destruction of tetrahedral local order by temperature, pressure or ionic solutes 
facilitates the coupling of local, librational modes to the network motions and therefore the 
parallel behaviour of the structural and diffusional anomalies is not surprising. The power 
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spectral analysis shows that the boundaries of the two regions cannot be identical because the 
structurally anomalous region is defined by the average values of the local order parameters 
whereas the diflusionally anomalous region is determined by the dynamical correlations in 
the fluctuations in the local order and energy. 

In this paper, we apply the power spectral analysis techniques described above to under- 
stand the diffusional anomaly in liquid silica, modeled using the BKS potentia l 23 ! 24 . While 
the behaviour of liquid silica and water has been shown to be parallel in many respects, 
there are some important difference o 8 i 10 i 25 i 26 i 27 i 28 . Among all the crystalline polymorphs of 
silica, only /3-cristobalite has a melting line with a slightly negative slope having a change 
in volume equal to -1% as compared to -8.9% seen for water—. From the chemical bonding 
perspective, silica is thought of as an ionic liquid with strong covalent character. The tetra- 
hedral motif in water is displayed by the Walrafen-pentamer structure, while in silica and 
many silicates, the SiC-4 tetrahedra are the crucial structural unit. The 0-H- • -O hydrogen 
bond in water is nearly linear while in silica, the Si-O-Si bond angle varies between 140- 
155°G 2 ^. It is therefore interesting to consider the extent to which the dynamical behaviour 
of the two networks is similar. Based on our studies of water, we have found the tagged 
particle (atom or molecule) energies as a convenient dynamical variable which is sensitive 
to local order and energy of the molecular environment. In addition, since it is defined as 
the interaction of a tagged atom or molecule with all other particles in the simulation cell, 
it is reasonably sensitive to network reorganisations over a fairly wide range of length and 
time scales. Therefore we focus on the tagged ion potential energies as the important local 
variable in a tetrahedral liquid and study the corresponding static distributions as well as 
the power spectrum of the fluctuations. 

The paper is organised as follows. Section 2 discusses the computation of the tagged ion 
potential energies when Ewald summation is used to evaluate the long range, Coulombic in- 
teractions. Section 3 summarises the aspects of power spectral analysis required in this work. 
Computational details regarding the potential model and molecular dynamics simulations 
is given in Section 4. Section 5 and 6 contain the results and conclusions respectively. 
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II. COMPUTATION OF TAGGED MOLECULE POTENTIAL ENERGIES 



Parametric potentials for silica are typically constructed as a sum of long-range, Coulom- 
bic interactions, U cou i, and short-range van der Waals interactions, U v d w , such that 

U = Ucoui + U vdw (4) 

The Coulombic contribution comes from the charge-carrying Si and O atoms. We assume 
that the short-range interactions include only two-body terms, but if three-body terms are 
present the computation can be adapted. The functional form of the short-range potential 
used in this work is discussed in section 4.1. For a pair additive system, if usi and uo are 
the tagged atom energies and ngi and no are the number of Si and O atoms respectively, 
then the corresponding ensemble averages are related to the average configurational energy 
as: 2(U)/N = (n St /N)(u St ) + (n /N)(u ). 

The short range van der Waals interactions contribute an amount (u^ dw ) to the tagged 
atom potential energy Ui. Since the interactions are short-ranged, the minimum image 
convention can be applied and u\ can be evaluated by summing over all pair interaction 
contributions between atom i and all other atoms j located in the central simulation cell. 
If the short-range potential has a finite value at the spherical cut-off radius, a density- 
dependent long-range correction term must be included in u v i dw . 

To evaluate the Coulombic contribution to the tagged molecule potential energy, one must 
consider the Ewald summation technique for evaluating long-range force a 30 ' 31 . The charge 
density distribution for a collection of N point charges located in the central simulation cell 
is given by: 

N 

P( r ) = H^( r - r *) (5) 

i 

Note that r is a 3-dimensional position vector and is the position of the i th charge and 
that the central simulation cell is set up to be electrically neutral. The simulation cell is 
taken as cubic with edge length L. The Coulombic interaction energy between all charges 
in the central simulation box and all other periodic images is given by: 

1 00 - - aa- 
[/ ml =TTV ^ r (6) 

where Yij is the vector distance between charges g; and qj, n = (n x L,n y L,n z L) denotes the 
set of possible translations of the simulation cell. The central simulation box corresponds 
to the n = case and the prime on the first summation indicates that the series does not 
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include the interaction % — j for n = 0. Unlike in the case of the short-range van der Waals 
interactions, the inverse distance dependence of the Coulombic interaction implies that the 
summation cannot be restricted to just the central simulation cell and the above series will 
converge very slowly. Ewald summation replaces the slowly convergent series in equation 
(5) by a sum of two rapidly convergent series obtained by replacing the charge distribution 
p(r) by a sum of two charge distributions, p screen and p ga uss- To obtain p SC reen(r), each point 
charge qi is surrounded by a neutralizing charge distribution, usually a Gaussian distribution 
centred at the location of qi. To ensure that p(r) = p S creen(?) + Pgauss^), Pgaussi?) must 
consist of Gaussian charges which neutralise the screening charges used to construct p SC reen- 
The potential energy contribution due to the screened charges is denoted by U screen . The 
advantage of choosing neutralising Gaussian distributions of the form: 

Pi( r ) = ex P(-^ 2 r 2 ) (7) 

is that U screen can be readily evaluated using the complementary error function, erfc, as 
follows: 

Uscreen 2^ 2^ 2^ A I- 1_ n | ^' 

z i=l j=l | n |=o * nt \'tj^ n \ 

This contribution is frequently referred to as 'real space' summation and its rate of conver- 
gence depends upon the parameter a of the neutralising Gaussians. The larger the value 
of a, the narrower the gaussian distribution and the shorter the range of interaction of the 
screened charges. Typically, a is chosen to ensure that the minimum image convention can 
be applied i.e. the range of the interactions is shorter than half the length of the simulation 
cell i.e. only the |n| = case needs to be considered. One can now rewrite: 

N N 

2 i=i r- 
1 N 

= ^E«r n (io) 

where u\ creen represents the contribution from the screened charges to the tagged atom 
Coulombic energy of atom i. 

The contribution of the Gaussian charge distribution, p ga ussij), to the Coulombic inter- 
action energy can be evaluated very simply in reciprocal or /c-space and is denoted by U rec . 
The electrostatic potential, Q ga uss, associated with this charge distribution can be evaluated 
in real space using: 

- V 2 $ sauss (r) = 4np gaU ss(r) (11) 



Uscreen ^ ] ^ ' . i i (9) 



qtqj erfc(a\ 
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or in the Fourier form in reciprocal space as: 

k 2 <5> gauss {k) =47cpg auss (k) (12) 

where p 9cm ss(k) is the Fourier transform of p(r). Fourier transforms of arrays of Gaussians 
can be computed analytically and one can then write 31 

U reC = \ Qi^gaussir) (13) 

Taking explicit summations over the partial charges, indexed by j, and the reciprocal lattice 
vectors, indexed by |n| one can write 

2 



1 ~ exp(-A; 2 /4« 2 ) 

Urec ~ 2^ 



2V e |n ^ P 



N 



^2qjexp(-ik ■ 



(14) 



j 

Using the above expression, U rec can be calculated using a double loop over |n| and j. 
Computing the tagged atom contribution, requires the modulus squared term to be replaced 
by a double sum over i and j, where i and j label the partial charges. Rearranging the above 
expression 

1 q- N 

Urec = TyTr YjYj exp(zfc • fj - k 2 / Aa 2 ) E q% exp(-ik • f n ) (15) 

1 q N 

= ^rEE T^exp(-A; 2 /4a 2 )^ftexp(-iA;(rj - rj) (16) 
z v o£ j fc ^o K i=o 

= ^E^ rec ( 17 ) 

Z 3 

Since the wj ec contribution multiplied by (—ik) gives the force acting on charged site j, the 
Ewald summation subroutines supplied with the DL_POLY program require only a small 
modification in order to extract ^ e c32 i 33 . 

It should be noted that the above expression for U rec includes a self-interaction term 
for the interaction of each Gaussian charge distribution with itself. When computing the 
Ewald summation to evaluate U cou i, it is necessary to remove these redundant self-interaction 
contributions, U se if, which are given by 

N „2 ( 

7T 



N a 2 a 
^/ = E< / = (V4 7 re)E% (18) 



3=1 3 

The complete expression for U cou i will be 



U C oui U rec -\- U screen U se ij (19) 



The contributions to the tagged potential energy that comes from the electrostatic interac- 
tions is given by [u\ ec + u\ CTeen — ). The overall tagged molecule potential energy, Ui, of 
molecule % will therefore be given by : 

m = uf w + < ec + uf reen - uf f (20) 



III. POWER SPECTRAL ANALYSIS 



The power spectral density of an observable A(t) as a function of time t, has been defined 
in equation (3). In this work, the dynamical observables are the tagged atom potential 
energies, Ui(t), where i corresponds to a Si or O atom. Since we sample the potential energy 
signal at discrete times, the definition of the power spectral density must be changed to 
allow for discrete sampling 20 : 



\F k 
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S k = S(f k ) = ^ (21) 

where N is the number of samples, and F k is the Discrete Fourier Transform (DFT) evaluated 
at the k-th frequency, i.e., 

N-l 

F k = F(f k ) = J2 A n e 2mnk / N (22) 

n=0 

(here we use the same normalization as ref. 20, the k-th frequency is f k = k/(NAt), At 
is the sampling interval, and k = 0, ...,N — 1). The dynamical time scales present in 
the system will clearly determine the behaviour of the power spectrum. Consider a simple 
situation in which there is a single relaxation frequency A in the system. If several such 
relaxation events take place at an average rate n at times t k , k = 1,2,..., then each such 
event will give rise to an exponential decay in some appropriate dynamical variable of the 
form A exp[—X(t — t^)] for t > tj, (and for t < t k ) leading to an overall time-dependent 
signal A(t) = J2k tk)- The power spectrum of the fluctuating part of this signal is^ 

s ^ - x^kjf < 23 > 

which is frequently referred to as Debye relaxation. In the limit / ^> A, the Debye term 
approaches a simple 1/ f 2 behaviour. 

If there are multiple time-scales, i.e. if the relaxation frequency of each of the single 
events is drawn from a probability distribution, g\(\), then the power spectrum must be 
computed as: 



If g\(\) corresponds to a uniform distribution over the range Ai to A2 and the amplitude of 
each pulse remains constant, then 

s{f) = x^xlL vT$^jf dX 

A ° U [arctan(A 2 /27r/) - arctan(A 1 /27r/)] (25) 



2tt/(A 2 -Ai) 

In the low-frequency limit such that f ^ X± ^ X2, the spectrum flattens and we observe 
white noise with (S(f)) = A^n. In the intermediate region where Ai <C / <C A2, we obtain 

s(/) = wS^T) (26) 

which corresponds to the l/f a regime with a — 1. In the high-frequency limit for which 
/ S> A 2 ^> Ai, we obtain a l/f 2 noise with S(f) = A\nj '(2ir f) 2 . Other distributions of 
time scales may also give rise to 1/ f a behaviour; for example, a distribution -P(A) oc A _/3 
over the same frequency interval will generate S(f) oc l/f 1+f3 . Ref.— provides a gallery of 
spectral densities associated with different models. In general, a spectral index a in the 
interval 0.5 < a < 1.8 may be taken as an indication of multiple time-scale dynamics while 
exponents very close to 2 indicate single time-scale behaviour. A separation of two decades 
in time scales between Ai and A 2 is sufficient to observe a 1//" regime with a close to 1. 
Unless otherwise stated, in this work, reference to l/f a behaviour will imply the multiple 
time-scale regime with a between 0.5 and 1.8, rather than the Debye tail with a ~ 2. 

In the case of a networked liquid, like water, the multiple time scale region spans approxi- 
mately two decades. There may be additional features, such as high-frequency resonances or 
low-frequency Debye relations. One can, in principle, disentangle the different contributions 
to the complex network dynamics using a model-based fitting procedure^. For the purpose 
of quantitatively characterising the changes in the hydrogen bond network dynamics in the 
diffusionally anomalous region, we have found it sufficient to identify the l/f a region from 
the In S(f) versus In / plot and obtain the scaling exponent a by numerical fitting. 

IV. COMPUTATIONAL DETAILS 
A. Potential Model 

We use the van Beest-Kramer-van Santen (BKS) potential to model interatomic interac- 
tions in silica since it has been extensively used to study liquid state anomalies and the glass 
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transition 2 ^!. The BKS potential is pair-additive and contains long-range Coulombic and 
short-range two-body contributions. The pair interaction between atoms i and j is given 
by: 

^BKsiTij) = -r^- + A ije xp- b ^ - ^ (27) 



1~ o/v, ' r. 



where r„ is the distance between atoms i and j carrying charges qi and qj and Ay, fry and 
Cij are the parameters associated with the Buckingham potential for short-range repulsion- 
dispersion interactions. In the original BKS potential, the Buckingham parameters for the 
Si-Si interactions are taken as zero. At high temperatures and pressures, the BKS potential 
exhibits an unphysical divergences in the interaction energy when the Si-0 distance becomes 
very small. To remove these divergences, a short range correction term was added by Poole 
et. al 26 to the BKS potential and the parameters adjusted in such a way that the original 
form was left unchanged at larger separations while the negative divergence was avoided at 
smaller separations. The modified pair interaction is 



Tij) = 4>bks + 4q 



\ 30 / 

?ii ] _ ( 111 



(28) 



where ey and cry are the energy and length scale parameters for the 30-6 Lennard- Jones 
interaction. The parameters for the modified BKS potential used in this work are given in 
Table I. 



B. Molecular Dynamics 

Molecular Dynamics simulations of a system of 150 Si and 300 O ions were carried out in 
canonical (N-V-T) ensemble, using the DL_POLY software packag o 32 ' 33 , under cubic periodic 
boundary conditions. The effects of electrostatic (long-range) interactions were accounted 
for by the Ewald summation metho d 30 ' 31 . The non-coulombic part was truncated and shifted 
at 7. 5 A. A Berendsen thermostat, with the time constant T£=200ps, was used to maintain 
the desired temperature for the production run. The Verlet algorithm with a time step of lfs 
was used to integrate the equations of motion. The system was simulated at 5 temperatures 
in the 4000K to 6000K range lying along 8 isochores in the density range from 1.8 g cm -3 to 
4.2 g cm -3 . For temperatures above 5000K, the system was equilibrated for 3-4ns followed 
by a production run of 5- 7ns. For lower temperatures, an equilibration period of 5-6ns 
was followed by a production run of 8- 10ns. The starting configuration for the liquid state 
simulations was generated by replicating the cubic unit cell of /3-cristobalite to generate a 
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cube of edge length 21.48v4 37 ' 38 which was then suitably truncated to obtain a simulation cell 
containing 150 Si and 300 oxygen atoms. The system thus obtained was heated gradually 
to form the equilibriated liquid state. Our simulation results are in agreement with previous 
result o 1Ql26 within the statistical error bars. 

Figure [1] shows the P(p) curve at different temperatures. It can be seen that at all den- 
sities below 3.0 g cm -3 , the pressures are negative for the temperatures studied, indicating 
that liquid silica is stretched rather than compressed. The isothermal compressibility, de- 
fined as kt = (1/ p)(dp/dP), is negative only for the p = 1.8 g cm -3 isochore. The spindonal 
density for which kt = is located around 2.0 g cm -3 . Below this density, the system is 
likely to be inhomogenous. In our simulations, the state points along the 1.8 g cm -3 iso- 
chore are thermodynamically unstable; however, to present a complete study to explore the 
diffusional anomaly we consider this density. 

C. Power spectral analysis 

In this study, we focus exclusively on the power spectra associated with tagged atom 
potential quantities, uo and Usi, which were stored at intervals of 10 fs during the MD 
runs. We represent power spectra of tagged potential energy fluctutations as S u (f) but 
while referring to S u (f) spectra of atoms of a specific type A, we use the notation Sa(/)- 
The sampling interval of 10 fs corresponds to a Nyquist frequency of 1666 cm -1 . The value 
of the Berendsen thermostat time constant (tb) provides the lower limit on the frequency 
range over which we can obtain reliable power spectra; thus, Tb = 200ps corresponds to a 
lower frequency limit of 0.165 cm -1 . Standard Fast Fourier Transform routines were used 
with a square sampling window^. The normalisation convention was chosen such that the 
integrated area under the S(f) curve equalled the mean square amplitude of the time signal. 
Windows containing 2 19 data points were used for Fourier transformation. Statistical noise 
in the power spectra was reduced by averaging over overlapping time signal windows as 
well as over individual tagged particle spectra. In a given frequency interval showing 1/ f a 
behaviour, linear least squares fitting of In S(f) was done to obtain the a values. 
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V. RESULTS 



Figure 2 shows the behaviour of the self-diffusivities, Dsi and Do of Si and atoms 
respectively, in BKS silica melt as a function of density along different isotherms. At any 
given state point, the oxygen self-diffusivities are greater than those of silicon atoms by 
approximately a factor of two but the density and temperature-dependent trends are very 
similar. Subsequently, in this paper, we will focus only on D$i since silicon atoms are the 
sites with the local tetrahedral symmetry. The diffusionally anomalous region is clearly 
demarcated in the isotherms lying between 4000K and 5000K. The density of minimum 
diffusivity, p m j n , is approximately 2.0 g cm -3 over this temperature range while the density 
of maximum diffusivity, p max , is approximately 3.4 g cm -3 . In the case of SPC/E water, 
the region of anomalous diffusivity is seen below 280K and the p m in and p ma x values are 
approximately 0.9 g cm -3 and 1.1 g cm -3 respectively. 

Figure 3 compares the power spectra generated by fluctuations in potential energies of 
tagged oxygen and silicon atoms in silica melt at 4000K along the 2.0 g cm -3 isochore with 
that of water at 0.9 g cm -3 at 230K. The chosen state point for both systems corresponds to 
the density minimum along the appropriate isotherm. Despite the difference in intermolec- 
ular interactions between a molecular liquid and an ionic melt, there is a striking parallel in 
the power spectral features suggesting a very similar pattern in the dynamical correlations. 

The Ssi{f) spectrum shows a sharp peak at 1100 cm' 1 corresponding to the to the Si- 
O-Si asymmetric stretch vibration seen in both the experimental as well as simulated IR 
spectru m 35 ! 36 . The power spectrum also shows a secondary, less intense peak at 800 cm -1 
corresponding to O-Si-0 bending. It should be noted that the tagged particle power spec- 
trum reflects the intrinsic dynamics of the liquid state network while the infra-red spectrum 
highlights those modes which couple to the electric field of the electromagnetic radiation in 
the IR frequency range. The localised high frequency peaks at 1100 and 800 cm' 1 correspond 
to well-defined local vibrational modes of the SiC>4 tetrahedra and mirror the librational peak 
in water which is, however, somewhat broader and more intense. 

The region between 0.2 and 400 cm -1 in Ssi{f) may be regarded as essentially a multiple 
time scale or l/f a region with a frequency- dependent value of a. At higher temperatures 
or densities, the very low frequency region shows a cross-over to white noise. The multiple 
time scale region of water spans a similar frequency range. 

The Ssi(f) and So(f) spectra are qualitatively very similar. The tagged oxygen power 
spectrum, So(f), typically has a lower intensity than the Ssi(f) power spectrum because 
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fluctuations in the tagged energy of the oxygen atoms are significantly lower than those for 
the Si atoms due to the relative magnitude of the charges. The high-frequency peaks in 
the So(f) spectra occur at essentially the same locations as in Ssi(f) but have significantly 
different intensities suggesting that tagged energies of the two types of atoms have different 
sensitivities to the local vibrational modes. 

Section 2 discusses the different contributions to the tagged atom potential energy i.e. 
van der Waals (u v d w ), screened charge (u screen ), self- interaction (u se if) and reciprocal space 
(u rec ). While u se if is a constant for a given atomic type, the power spectra associated with 
fluctuations in the other three quantities at the state point (T=4000K, 2.0 g cm -3 ) is shown 
in Figure 4. The behaviour of S u (f) is almost exactly reproduced over much of the frequency 
range by the screened charge contribution. The spectral power associated with the van der 
Waals and reciprocal space contributions is much lower and therefore has a minor effect on 
the overall spectrum except in the high-frequency region. Unlike in the case of water, the 
u v dw contribution produces a relatively unstructured power spectrum, specially in the case 
of oxygen, while the Coulombic interactions are responsible for much of the high-frequency, 
short length scale structure. This is corroborated by the static distributions of the various 
contributions to the tagged particle potential energies shown in Figure 5. The width of the 
P(u rec ) and P{u v dw) distributions is much narrower than that of P(u screen ) both in the case 
of Si and O atoms. 

To understand the effect of pressure and temperature on the tetrahedral network in 
liquid silica, we examined the power spectra as well as static distributions of tagged particle 
potential energies as a function of temperature and pressure. Figures 6 and 7 show the 
results along the 2.0 g cm^ 3 isochore and the 4000K isotherm. 

Figure 6(a) shows that, along a given isochore, the crossover to white noise occurs at 
increasingly higher frequencies as temperature increases. Moreover the local vibrational 
peak broadens and the upper bound of the frequency range for multiple time-scale or l/f a 
behaviour increases. This progressive coupling of local modes into network vibrations leads 
to a change in the exponent a of the 1/ f a region. At the lowest temperatures, the frequency 
dependence of a is more pronounced. For example, at 4000K and 2.0 g cm -3 one can identify 
a high frequency regime from 20-200 cm -1 with a = 0.8029 while between 1-20 cm -1 , the a 
value is 1.6526. Since we wish to use a as a measure of the degree of coupling of local and 
network modes, we consider a in the 20-200 cm -1 region. As we show below, this change is 
quantitatively correlated with increasing frequency in the anomalous region. 
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The effect of increasing density or pressure on the power spectrum is shown in Figure 
6(b). The crossover to white noise occurs at about 1 cm -1 for densities of 3.0 g cm -3 or 
more. The effect of increasing pressure results in loss of the bimodal character of the high- 
frequency peak. The peak at approximately 800 cm -1 corresponding to the O-Si-0 bend 
is clearly much more sensitive to pressure than the Si-O-Si asymmetric stretch; this is not 
surprising since it is known that compaction results in distortion of bond angles resulting 
in more efficient packing of SiC>4 tetrahedra. This is different from the loss of the distinct 
identity of the librational peak seen in water with increasing density or temperature. 

Figure 7 shows the static distributions, P(u) of the tagged atom potential energies. At 
low temperatures or low densities, the P(us%) distributions are multimodal with six dis- 
tinct peaks. This suggests that the energetic environment of a given Si atom depends on 
bond lengths and bond angles associated with the nearest four oxygen atoms forming the 
local Si0 4 tetrahedral unit. When all the angles and bond lengths are optimal, the lowest 
energy environment is obtained. The presence of distinct peaks suggests that there is a 
small energetic barrier to transit from one energetic environment to another. The effect of 
temperature is to allow for more facile transitions between different energetic environments 
which results in a loss of multimodality but no additional broadening of the distributions. 
In contrast, the effect of pressure is to both broaden as well as destructure the distributions. 
The P{uo) has two major and one minor peak at low temperatures or low densities. With 
increasing temperature, the minor peak becomes part of the shoulder of the low energy peak. 
With increasing density, the minor, central peak disappears and the overall wiidth of the 
distribution broadens. 

The earlier study on water has established the strong correlation between the diffusivity 
and the scaling exponent a of the multiple time-scale regime in the anomalous regime. 
Figures 8 and 9 show that this relationship also holds in the case of liquid silica. Figure 8 
shows the behaviour of the scaling exponent a^ 1 of the Ssi(f) spectra as a function of density 
along different isotherms. The parallel with the diffusivity plot in Figure 2 is obvious. It 
should be noted, however, that a is essentially constant at 1.5 at high temperatures and 
densities when the liquid is in the diffusionally normal regime. Figure 9 shows a correlation 
plot of af l with In D si where the clustering of a values around 1.5 can be seen. 
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VI. DISCUSSION AND CONCLUSIONS 



This study demonstrates the strong parallels in the dynamical behaviour of the liquid 
state networks in silica and water in the anomalous regime. Power spectral analysis is a 
convenient tool for characterising the complex dynamics of network-forming liquids, such 
as silica, since the multiple time-scale behaviour of the network is captured by the range 
and exponent of the l/f a regime. It is particularly useful for understanding diffusion in 
the random network of the liquid. In normal liquids, diffusion proceeds via collisions and 
relaxations of the local nearest-neighbour cage. In networked liquids, local vibrational modes 
can couple to network reorganisations, thereby facilitating diffusion. The details of the 
process in water and silica differ, but in both cases, the degree of coupling can be indexed 
by the exponent a of the multiple time-scale regime. As a consequence, the diffusivity is 
strongly correlated with the scaling exponent, specially in the anomalous regime. As the 
network connectivity is attenuated by temperature and pressure in the normal regime, the 
correlation begins to break down. 

The results of this study provide some insights into the microscopic connection between 
the structural and diffusional anomalies. The boundaries of the structurally anomalous 
region are determined by extrema in the local order parameters. The progressive destruction 
of local tetrahedral order by temperature or pressure has an effect on the local atomic or 
molecular environment, as can be seen from the static distribution of the tagged particle 
potential energies. As the local tetrahedral order is destroyed, the vibrational coupling 
of local modes, specially the O-Si-0 bend, to the network reorganisations becomes much 
stronger causing a rise in diffusivity. The diffusivity does not, however, rise indefinitely with 
increasing compression because for sufficiently large values of the density, the extent and 
connectivity of the network is severely attenuated and the liquid shifts from the anomalous to 
the normal regime. Therefore one would expect the diffusional and structural anomalies to 
be strongly correlated. Whether the structural anomalies precede or succeed the diffusional 
anomaly (as in water and silica respectively) must depend on the details of the potential 
energy surface, such as the nature and extent of vibrational coupling. Based on the similarity 
in power spectral behaviour of SPC/E and TIP5P-E models for water 19 , we would expect 
different effective pair potentials for silica to have very similar qualitative behaviour. It is 
possible, however, that three-body potentials for silica will show more significant deviations, 
such as a reversal in the relative order of the two types of anomalies 35 . 

The parallel behaviour of power spectra of water and silica in the diffusionally anomalous 
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region suggests that all tetrahedral liquids should display this generic pattern of behaviour. 
It would be interesting to perform a power spectral analysis of other liquids showing water- 
like anomalies, such as two-scale ramp potentials 39 , to see if the detailed dynamical behaviour 
resembles that of tetrahedral liquids. 
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Table I 

BKS (modified) potential parameter o 23 ' 24 . 



Aij bij ^ij ^ij ^ij 

(kJ mol" 1 ) (A" 1 ) (kJ mok 1 ) (kJ mok 1 ) (A) 



- 134015 2.76 16887.3 0.101425 1.7792 
Si-0 1737340 4.87 12886.3 0.298949 1.3136 
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Figure Captions 

1. Equation of state for BKS silica: Pressure (P) as a function of density (p) along 
different isotherms. Systems at a density of 1.8 g cm" 3 are likely to be inhomogeneous. 

2. Dependence on density, p, of (a) self-diffusivity of Si, D S i, and (b) self-diffusivity of 
O, D along different isotherms. 

3. Comparison of power spectra associated with temporal fluctuations in tagged particle 
potential energy of Si abd O atoms at 4000K and 2.0 g cm" 3 with water at 230K and 
0.9 g cm" 3 . The power spectrum of O atoms has been scaled to half for clarity. 

4. Contributions of reciprocal space and screening from Ewald summation and van der 
Waals interactions to the power spectra of total tagged potential energy fluctuations 
of (a) Si at 4000K and 2.0 g cm" 3 and (b) O at 4000K and 2.0 g cm" 3 . Relative 
magnitudes kept unaltered. 

5. Contributions of reciprocal space and screening from ewald summation and van der 
Waals interactions to the static distribution of tagged particle potential energy of (a) 
Si at 4000K and 2.0 g cm" 3 and (b) O at 4000K and 2.0 g cm" 3 . In part (a) van der 
Waals and reciprocal space contributions, P(u v d w ) and P(w rec ) respectively, have been 
scaled by a factor of 1/10. In part (b) P(u vdw ) has been scaled by 1/4 while P(w rec ) 
has been scaled by 1/10. 

6. Power spectra of tagged particle potential energy fluctuations, S u (f), of Si in liquid 
silica along (a) 2.0 g cm -3 isochore from 4000K to 6000K in steps of 500K and (b) 
4000K isotherm at 1.8, 2.0, 2.2, 2.6, 3.0, 3.4, 3.8 amd 4.2 g cm" 3 . Relative magnitudes 
of power spectra for different temperatures and densities were adjusted (scaled) for 
clarity. 

7. Static distributions of tagged particle potential energies (in kJ mol" 1 ) of (a) Si along 
2.0 g cm" 3 isochore, (b) Si along 4000K isotherm, (c) O along 2.0 g cm" 3 isochore and 
(d) O along 4000K isotherm. 

8. Multiple time-scale exponent, function of density along different isotherms. 
The exponent af l is evaluated using the high frequency 1/ f a regime of 20 - 200 cm -1 
at 4000K and 4500K and from the point of crossover to 200 cm" 1 for 5000K, 5500K 
and 6000K. 
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9. Correlation plot showing In D S i against a„ for different isochores. The frequency 
range for evaluating af* are the same as in Figure 8. Units of diffusivity are taken as 
cm 2 s _1 . 
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